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ABSTRACT 

Many extrasolar planetary systems containing multiple super-Earths have 
been discovered. N-body simulations taking into account standard type-I plane- 
tary migration suggest that protoplanets are captured into mean-motion resonant 
orbits near the inner disk edge at which the migration is halted. Previous N-body 
simulations suggested that orbital stability of the resonant systems depends on 
number of the captured planets. In the unstable case, through close scattering 
and merging between planets, non-resonant multiple systems are finally formed. 
In this paper, we investigate the critical number of the resonantly trapped plan- 
ets beyond which orbital instability occurs after disk gas depletion. We find that 
when the total number of planets (N) is larger than the critical number (iV CI . it ), 
crossing time that is a timescale of initiation of the orbital instability is similar to 
non-resonant cases, while the orbital instability never occurs within the orbital 
calculation time (10 8 Kepler time) for N < N crit . Thus, the transition of crossing 
time across the critical number is drastic. When all the planets are trapped in 
7:6 resonance of adjacent pairs, N crit = 4. We examine the dependence of the 
critical number of 4:3, 6:5 and 8:7 resonance by changing the orbital separation 
in mutual Hill radii and planetary mass. The critical number increases with in- 
creasing the orbital separation in mutual Hill radii with fixed planetary mass and 
increases with increasing planetary mass with fixed the orbital separation in mu- 
tual Hill radii. We also calculate the case of a system which is not composed of 
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the same resonance. The sharp transition of the stability can be responsible for 
the diversity of multiple super-Earths (non-resonant or resonant), that is being 
revealed by Kepler mission. 

Subject headings: Celestial mechanics; Planetary dynamics; Planetary formation; 
Resonances, orbital 
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1. Introduction 

Seventeen years searchers of extrasolar planets have found more than 100 super-Earths 
and hot Neptunes (< 30Af e ). Some of them form multiple planet systems. Kepler 
Mission reports 885 of those multiples in 361 systems (Batalha et al. 2012; Borucki et al. 
2011a, 2011b). Period ratios of pairs of planets show some peaks at commensurable ratios 
(Lissauer et al. 2011; Fabrycky et al. 2012). Since orbital angles such as an argument of 
pericenter are unknown in many of the systems, it is not clear that the planets are actually 
in the mean-motion resonances. In data of first four months of Kepler Mission, Among 
158 multi-super Earth candidates and multi-Neptune candidates, about 25% planets have 
orbital periods which are almost commensurable with neighboring planets within 3% period 
ratio. Veras and Ford (2012) identified 70 non-resonant KOI (Kepler Objects of Interest) 
near resonant pairs which are not in a mean-motion resonance. Terquem and Papaloizou 
(2007) proposed a mechanism to form resonant planets near the disk inner edge. Possible 
mechanisms for small derivations from the resonance were also raised by Papaloizou and 
Terquem (2010). 

When planets are in the mean-motion resonance, the conjunction periods of the planets 
are expressed as the integer ratio of periods of each planet. That means, conjunctions of 
the planets always occur the same relative positions. Planets in a resonance can become 
stable depending on the configuration of conjunctions. For example, Neptune-Pluto are in 
3:2 mean-motion resonance. Although the orbits of them cross, they always avoid close 
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approach. Neptune-Pluto system are long-term stable (Cohen and Hubbard 1965). 

In a gas disk, growing protoplanets migrate toward their central stars due to type-I 
migration. The migration is stopped when a protoplanet arrives at the inner edge of the gas 
disk, often assumed to be at the corotation point of the star (~ 0.05 — 0.1 AU), if it exists. 
N-body simulations (Terquem and Papaloizou 2007; Ogihara and Ida 2009) showed that 
in the case of standard type-I migration (e.g., Tanaka et al. 2002) subsequently migrating 
protoplanets are usually not trapped in a mean-motion resonance at a first encounter 
with the protoplanet at the disk edge. After some close scatterings and collisions, they 
are eventually trapped in resonances because they are subject to relatively slow type-I 
migrations near the edge and eccentricity damping after relaxation. Through merging of 
many planets that have migrated to the inner edge, only several merged bodies finally 
remain in mean-motion resonant orbits. Although several inner planets are pushed inside 
of gas disk edge (< 0.05 AU) and others are in the gas disk (> 0.05 AU), they keep the 
relation of mean-motion resonances. These planets are spaced by 5 — 9 Hill radius with 
each other and stay stable even after gas depletion in which eccentricity damping no more 
operates. 

However, Ogihara and Ida (2009) have found that in the case of migration that is 
slow compared to the rates predicted for standard type-I migration, orbital evolution is 
totally different and final orbital configuration is non- resonant. In this case, subsequently 
migrating protoplanets are trapped in mean- motion resonances. As a result, about 40 
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small protoplanets queue in low order mean-motion resonances having closer separations 
in the gas disk at > 0.05 AU. Few of inner planets are pushed into the inner cavity, 
because of "eccentricity trapping" caused by torque balance between migration torque and 
edge torque (Ogihara et al., 2010). In contrast to the fast migration case, the planets 
become orbitally unstable after the gas depletion, i.e., their eccentricities are excited and 
their orbits start crossing. Through collisions and merges of planets, they are kicked out 
of the resonances and finally several planets are formed in non-resonant orbits with the 
large orbital separation (~ 15 — 20 Hill radius). Although resonant systems are generally 
more stable than non-resonant systems, results of Ogihara and Ida (2009) showed that 
multi-planet systems whose planets are initially in overpopulated resonant orbits become 
unstable in relatively short timescale after gas depletion and end up with dynamically 
relaxed non-resonant systems. 

Although the crossing time at which multi-planet systems have been extensively studied 
in gas free environment by N-body simulations (e.g., Chambers et al., 1996; Yoshinaga 
et al., 1999; Zhou et al., 2007), they only investigated non-resonant systems. Because 
of type-I migration and eccentricity damping, the systems that we consider are deep in 
resonances, so that the crossing time can be very different from that found by the previous 
studies. 

Although the previous studies on crossing times were concentrated on non-resonant 
systems, the observed resonant multi-super Earths' systems near stars and N-body 
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simulations suggest the occurrence of resonance trapping as a consequence of planetary 
migrations. In this paper, we mainly consider high-integer resonance e.g., 6:5 or 7:6. 
This is because previous N-body simulations suggest that proto-super-Earths are once 
in a close resonances with separations of ~ 5 — 9 Hill radius and cause instability due 
to overpopulation later on. We calculate the crossing time of systems in first-order 
mean-motion resonances, by changing the total number of plants, the orbital separation in 
mutual Hill radii, and planetary masses. In §2J we summarize previous studies of N-body 
simulations to evaluate the crossing time, t croS s- In 331 we explain numerical models of our 
simulations. As commensurability of orbital periods does not necessarily mean mean-motion 
resonance, in §U we study both cases that resonant and non-resonant systems having the 
same orbital commensurability. We discuss the results in §5j 



2. Previous Studies on Crossing Time of Multi-Planet Systems 

Here we summarize previous studies of N-body simulations to evaluate t cross for 
non-resonant systems in order to make clear the purpose of our simulations. Chambers 
et al. (1996) first investigated the crossing time of multi-planet systems, at which the first 
close encounter occurs. They performed orbital calculations of equal mass protoplanets 
with various mass from 1O _9 M to 10 _5 M Q . They put planets on initially circular and 
coplanar orbits with mutual separations aj — a, = Krmj, where r-mj is the mutual Hill 
radius of planets i and j, setting the innermost planet at a± — 1 AU. They repeated orbital 



-8 - 

simulations 3 times for the same Hill separation K changing planetary longitudes. The 
calculations were continued until the first encounter, which is defined by the time when 
distance between two planets becomes smaller than one mutual Hill radius occurs. They 
found that t CT oss is given approximately by an empirical relation, 

log Across = bK + C, (1) 

where b and c are constants. When systems are composed of 3 planets whose mass are 
lO~ 7 Af , b = 1.176 and c = 1.663 for example. The values of these constants depend on 
planetary mass (M p ) and the number of planets (N). But when a system has more than 5 
planets, adding further planets does not make significant difference to the stability of the 
system. 

The crossing time of protoplanets also depends on initial eccentricities and inclinations 
of protoplanets (Yoshinaga et al. 1999; Zhou et al. 2007). Yoshinaga et al. (1999) found 
that the two constants b and c decrease proportional to root mean square of eccentricities 
and inclinations. The dependence of constants b and c on planetary mass, M p , was studied 
by Duncan and Lissauer (1997) and Zhou et al. (2007). Duncan and Lissauer (1997) 
studied the crossing time of Uranian satellite system with multiplied satellite mass and 
Zhou et al. (2007) investigated the crossing time of protoplanetary systems whose settings 
are similar to Chambers et al. (1996). From numerical calculations of the crossing time 
with different masses, they concluded logt croS s oc logM p . They also empirically expressed 
dependence on eccentricity. The crossing time of systems containing retrograde planets 
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were studied by Smith and Lissauer (2009). They have found that systems with mixture of 
planets in retrograde and prograde orbits are more stable than the system which has only 
prograde planets provided that the total number of planets and their orbital separations in 
mutual Hill radii are the same. 

When planets embedded in protoplanetary gas disk, drag force which damps 
eccentricities also affects crossing time (Iwasaki et al., 2001, 2002; Iwasaki and Ohtsuki 
2006). When the crossing time without drag force is shorter than eccentricity damping 
timescale, the drag force hardly changes the crossing time. Otherwise, the orbital crossing 
time is at least 200-times longer than that without the drag. This result implies that 
multiple planet systems do not start orbital instability until disk gas is sufficiently depleted. 
If the effects of type-I migration and disk inner edge are taken into account, the systems 
can become resonant. We will show that the resonant configuration stabilizes the systems 
even after disk gas is completely depleted, if the number of planets is smaller than a critical 
value. 



3. Numerical Model 

We consider a situation that planets are brought to current locations near their 
host star by type-I migration, which leads the systems to resonant configuration in the 
protoplanetary disk. Planetary growth simulations including type-I migration and disk 
inner edge by Ogihara and Ida (2009) suggest that similar-sized planets are trapped in 
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the resonances. We consider that planets have equal masses (M p = 3M e — 30M e ) and 
coplanar orbits around the central star with M* = 1M Q in all cases (non-zero inclination 
cases are studied in Yoshinaga et al. 1999). Using 4th-order Hermite scheme, we continue 
calculations until a distance between planets becomes smaller than their mutual Hill radius 
for at least one pair or until we reach to an upper limit of orbital evolution time. We set 
the upper limit of our calculation at 10 8 Kepler time of the innermost planet (a\ = 0.1 AU 
for all cases). Planetary mass (M p ), orbital separation normalized by mutual Hill radius 
(K), and the total number of bodies (N) are treated as parameters. 

In this paper, we target on the first-order mean- motion resonances, i.e., planets have 
p + l:p period relation, according to the results by Ogihara and Ida (2009). Planets which 
are in a mean-motion resonance have a relation between their pericenters and a point of 
conjunction. Even if whether planets have a periods ratio of p + q:p, it does not guarantee 
that they are in a mean-motion resonance. Thus, for the same Hill separations and 
planetary mass, we can set up both resonant and non-resonant systems. In the non-resonant 
cases, planets have the p + l:p period ratio, but their initial longitudes are given randomly. 
In the resonant cases, we put planets in the mean-motion resonance using orbital migration. 
The orbital migration automatically leads the planets to the mean-motion resonance with 
appropriate resonant angles. To extract the effect of a mean-motion resonance on the 
orbital stability, we compare the crossing time for the resonant cases with that for the 
non-resonant cases. 



3.1. Configuration of System 

Here, we explain how we control separations of neighboring planets in the same 
resonance using only one parameter K. Mutual Hill radius of the z-th planet and the 
(z + l)-th planet is given by 

fMi + Mi+A 1/3 (M % a t + M i+1 a i+l \ 
rHl ' l+1 ~ { 3Af, ) { M l + M l+1 ) ' (2) 

where a, is z-th planetary semi-major axis, Mi is z-th planetary mass, and M* is stellar 
mass. Using a factor K, the orbital separation of neighboring planets is expressed as 

Oi+i ~ a i = Kr m , i+1 . (3) 

In our simulations, all planets have an equal mass Mj = M i+ i = M p and neighboring 
plantes have period ratio of p + q:p, i.e., a i+ i/di = (p/p + q)~ 2 ^ 3 . Using these relations, K is 
expressed as 

2{(p + q)^-p^} (3MA 1/3 

( p + g )2/3 +p 2/3 \2M V ) ' 1 ' 

Note that K in the same for all the adjacent pairs for given p and q. The values of K 
we use are shown in Table [TJ Semi-major axis of the z-th planet (z > 2) in the first-order 
mean-motion resonance is expressed as 



2M-K 

where M = (2M p /3M,)" 1/3 . We use a x = 0.1 AU in all cases 



. 2M + K . 

a-i = 7tt> — ai, (5) 



3.2. 



Non-Resonant Cases 



For non-resonant cases, we set planets according to Eq. [5] without any special treatment 
like resonant case ( §3.3p . We integrate the equations of motion, 



where i refers to the i-th protoplanet (z = 1,2,---, N), G is the gravitational constant, and 
Vij is relative distance of the planet % and j. We perform 100 runs for each value of N. 



We form exact resonance situations by orbital migration simulations in a gaseous disk. 
After all planets are locked in a resonance, we gradually deplete the gas. We calculate 
crossing times of 5 resonances. Choices of the resonance and planetary mass are in Table 
1. By the choice, the Hill separation K is automatically adjusted as we explained in §3.11 
In easel and case2, we use typical values of resonant planets obtained in simulations by 
Terquem and Papaloizou (2007) and Ogihara and Ida (2009). The other sets are chosen to 
study how crossing time changes by these parameters. We calculate 10 runs in 6:5 and 7:6 
mean-motion resonance, and 3 runs in other cases slightly changing initial longitude. We 
follow Ogihara and Ida (2009)'s settings. The basic equations to form initial conditions are 





3.3. 



Resonant Cases 



dt 2 



GM^ - GM p ^f - £ GM p ^f + F damp + F 



(7) 
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where -Fdamp and -F m i g are the specific forces owing to eccentricity damping due to tides 
from the gas disk as a drag force (e.g., Tanaka and Ward 2004) and type-I migration, 
respectively. These force are given by Ogihara and Ida (2009) as 

F — = (t)(?) 4 K) fi[{0 ' 114K " rfi)+0 ' 176 " r}e ' 

+{-1.736(u fl - rfi) + 0.325v r }e e + {-1.088v, - 0.871zfi}e z ] , (8) 
F„ g = -ilT/.gg)'^.., (9) 

where f m is a scale parameter corresponding to an uncertainty in type-I migration. These 
additional forces arise from interaction with the gas disk. For the disk model, we use 

Sg = 2400 Miau) gcm ~ 2 ' (10) 

, r v -i/4 / r \ 1 / 8 

c - = 10x1 ° 5 (iau) [rj cms " 

where S g is a surface density of the gas disk and c s is the sound velocity. During migration 
simulations, the surface density is 1.4 times larger than that of the Minimum Mass Solar 
Nebula model, i.e., f g = 1. The sound velocity is that for an optically thin disk. We assume 
that the disk surface density smoothly vanishes with a hyperbolic tangent function at inner 
edge as 

where r e( j ge is a heliocentric distance of the inner edge of the gas disk, and Ar represents 
typical width of the inner edge. We choose r e d ge = 0.1 AU and Ar = 0.001 AU. In our 
calculations, we put all planets slightly outside of p+ hp resonance. Planets slowly migrate 
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inward and are trapped at p + l:p resonance automatically. After planets are captured in 
mean-motion resonances, we perform orbital integration for crossing time, decreasing gas 
density. For adiabatic gas depletion, we decrease / g as 

/g = «*p(-— \ (13) 

where Td ep means depletion timescale. The timing t = is the starting time of gas 
depletion. For adiabatic gas depletion, we normally take Td cp = 10 4 yr following Ogihara 
and Ida (2009). We check the dependence of r^ cp on the crossing time by changing it to 
r dcp = 10 4 yr, 10 3 yr, 10 2 yr in §4.1.31 The system is stable over t drag (defined in Appendix 
[A]) due to eccentricity damping. To distinguish the resonant effect from the stabilization 
effect due to eccentricity damping, we choose these timescales as Td cp . 

The timescale of eccentricity damping obtained by linear calculation of tides from the 
gas disk (Tanaka and Ward 2004) is 

^~^(i&Y^)W*(£T" (14) 

where M and L Q are solar mass and luminosity. The timescale of semi-major damping by 
standard type-I migration (Tanaka et al. 2002) is 

'---^(^)^^) 3/2 (3 V2 (fe) 1/4 - < i5) 

To stop the innermost planet at the disk edge of a\ = 0.1 AU and to trap following planets 
in a resonances, we adopt slow enough migration, / m > 50 (Ogihara et al. 2010). When a 
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planet is in the first-order resonance, the timescale of resonant libration is 



'lib 



P 



aF D 



-3.94613 



11.7 



xEU /sin 2 (<p /2)) yr, 



-1/2 



io- 



-1/2 



a \3/2 / M p 
0.1ALV I 10~ 5 M« 



-1/2 



d 

-2p - a— 

da 



U l/2> 



(16) 
(17) 



where a is the ratio of semi-major axis of the inner planet of the resonant pair divided by 
semi-major axis of the outer planet, is Laplace coefficient, E is the complete elliptic 
integral of the first kind, and (fio is the maximum width of resonant libration (Murray and 
Dermott 1999). For adiabatic gas depletion, gas depletion timescale Td ep should be much 
longer than t^. The gas depletion in this work is always adiabatic. 



4. Results 

4.1. Dependance of Crossing Time on Number of Planets 

4-1.1. typical orbital evolution 

In resonant cases, we make initial conditions by damping of e and a, as explained in 
§3.31 Fig. [T] shows evolution of the semi-major axes of planets. The planets migrate inward 
due to type-I migration and are trapped in 7:6 mean-motion resonances at t ~ —4500 yr. 
Although planets are still subject to type-I migration, the innermost planet is caught in the 
disk edge at 0.1 AU and the other planets do not migrate furthermore. After gas depletion, 
planets start instability at t ~ 73000 yr. We confirmed the mean-motion resonance from 
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plain commensurability by resonant angle. When a planet is in 7:6 mean- motion resonance, 
the resonant angle of outer planet is written as ip[ = 7Aj + i — 6 A; — z&i+i, where A is the 
mean longitude and w is the longitude of pericenter and the resonant angle of inner planet 
is written as tpi = 7Aj + i — 6Aj — Wi. In the case of systems that are formed by migration 
(Fig. [2]), resonant angles of inner planet librate around ip = and resonant angles of outer 
planet librate around <p' = tt. Since tp — and tp' = tt, all conjunctions occur when inner 
planet is at its pericenter and outer planet is at its apocenter. This means that only one 
configuration is possible when tp = and ip' = it. The resonant angles start circulation 
at about 70000 yr, just before the occurrence of instability. On the other hand, if we just 
put planets at semi-major axes of 7:6 resonance without such special treatment migration, 
resonant angles circulate and do not take a particular value (Fig. [3]). That is, this initial 
condition is non- resonant. 



4-1.2. non-resonant case 

In the non-resonant case, we put planets (M p = 1O _5 M ) on the orbital sep- 
aration of Aa = Oj + i — dj = 6.450rHi,i+i (K = 6.450). We calculate iV = 
3, 5, 8, 9, 10, 11, 20, 30, and 50 cases. Changing initial longitude randomly, we 
repeat the simulation 100 times for each N. Fig. |4] shows crossing time of 6:5 case (easel) 
having the different number of planets. The crossing time is normalized by Kepler time 
of the innermost body at a± = 0.1 AU. The circles are the crossing time of non-resonant 
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system. The solid curve is a least square exponential fit for the results of non-resonant 
system. Although there is some fluctuation, the crossing time decreases with N. However, 
in the region of N > 10, crossing time is almost constant, t cross ~ 10 4 Xk ep - This tendency is 
consistent with the result of Chambers et al. (1996) using 10 _7 M Q planets. We formulate 
the crossing time as a function of the total number of planets such as 

logtcross = gexp (-log AT) + h = + h, (18) 

where g and h are constants. From a least square fit for the data in Fig. HJ g = 6.21 and 
h = 3.39 (the solid curve of Fig. H]). We show in §4.21 how these constants depend on the 
Hill separation and planetary mass. 

4-1.3. resonant case 

The results in resonant case (easel) are summarized in Fig. |H Triangles, squares, 
and crosses are the crossing time with Td cp = 10 4 yr, 10 3 yr, and 10 2 yr, respectively. 
Three dotted curves are gas stabilization timescale (tdrag) which is defined in Eq. IA1I for 
r dcp — 10 4 , 10 3 , and 10 2 yr from the top. Symbols shown on the top horizontal axis 
indicate lower limits of crossing time, since crossing has not been detected within 10 8 Xkep- 
In Fig. HJ the crossing time of the resonant case is longer than that of the non-resonant 
case. Since the resonant cases are stabilized by resonant effect and eccentricity damping by 
gas, we have to estimate the timescale of the stabilization by gas to understand resonant 
effect. Gas drag is needed to form resonant systems, but its removal is required for the 
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evaluation of the crossing time. Since rapid gas depletion makes system unstable, adiabatic 
gas depletion is needed. Due to the adiabatic depletion, planets are stabilized by gas on 
timescales ~ tdrag (defined in Appendix [A]) . In Fig. [4j t cross ~ tdrag for N > 10. The gas 
drag effect is that t C ross cannot be shorter than tdrag suggested by Iwasaki et al. (2002). The 
crossing time in large N is dependent on Td cp . As Td ep becomes shorter, we can diminish the 
gas drag effect and t cr0 ss approaches that in the non-resonant case for N > 10. Then, we 
find that t cross jumps up by several orders of magnitude at iV ~ 8. This jump-up is due to 
resonant effect. In the following, we examine the critical total number of planets at which 
t cross increase abruptly with decreasing N . Note that iV cr it is almost independent of Td cp . 
Since resonant libration timescale is about 20 yr in easel, small Td ep is not long enough to 
guarantee adiabatic gas depletion in these cases. Therefore, we choose Td cp = 10 4 yr in the 
following, although an off-set of ~ idmg for Td ep = 10 4 yr is added in the results. 

4.2. Orbital Separation and Mass Dependence 

In this subsection, we show dependence of crossing time on the orbital separation in 
mutual Hill radii and the planetary mass. First, we change the orbital separation factor K 
fixing planetary mass. It corresponds to a change of p of the resonance. The results for 6:5 
resonance (K = 6.450), 7:6 resonance (K = 5.456), and 8:7 (K = 4.727) are shown in Fig. 
HI O and respectively. In non-resonant cases, g and h of Eq. [TS] increase with increasing 
K. This tendency is consistent with Chambers et al. (1996). The crossing time of resonant 
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planets shows discontinuity at a certain value of N as explained in §4.11 for 6:5 resonance 
case. The critical number of holding planet is N cr - lt = 4 for 7:6 resonances and iV cr i t = 3 for 
8:7 resonances. We find a tendency that N CT i t decreases with increasing p value of p + l:p 
resonance, i.e., decreasing the orbital separation in mutual Hill radii, as we show it in Fig. 

LH 

Next, we change planetary mass with fixed the orbital separation in mutual Hill 
radii. Since mutual Hill radius depends on planetary mass, we can choose some resonances 
having similar K by changing planetary mass. We choose planets of mass 10 _4 M o for 
4:3 resonances in case4. In this case, the orbital separation in mutual Hill radii is equal 
to K — 4.715 which is nearly equal to the case3 (8:7 resonance with M p = 1O" 5 M ). 
These results are plotted in Fig. |HJ In non-resonant cases, crossing time is longer for 
larger masses. According to Chambers et al. (1996) and Zhou et al. (2007), the crossing 
time of non-resonant systems increases with increasing planetary mass for K > 4 as long 
as planetary eccentricities are small enough. This tendency is also found in our resonant 
results. In resonant cases, iV cr i t is 7 in case4. Since iV cr it = 3 in case3, A^ crit decreases with 
increasing p value of p + l:p or decreasing planetary mass. 

Quillen (2011) suggests that three-body resonance overlap affects crossing times of 
non-resonant systems. It would also be the case in our resonant systems. Our results show 
that iV CI .i t increases with decreasing p value. It is consistent with the fact that the resonance 
overlap less occurs with small p value. 



-20 - 

We find that a pair which causes instability is always the nearest neighbors in outer 
region for resonant cases. It is related to our exponential depletion of gas density expressed 
in Eq. [13j Outer planets can cause instability while inner plants are still stabilized by gas. 
Planetary pair which causes instability depends on gas profile and the way of gas depletion. 

4.3. Chain of Resonance 

We have been studying crossing time of planets which are in the same p + l:p 
resonance for all the adjacent pairs. But there is not the cases for resonant planets that 
are observed and formed by N-body simulations of Ogihara and Ida (2009). For example, 
4 planets observed in KOI-730 have periods of 7.38469, 9.84978, 14.7845, and 19.72175 
days, respectively (Lissauer et al. 2011). Lissauer et al. (2011) suggests that these 
planets are in a chain of resonance of 8:6:4:3. In this system, planets are in the first-order 
mean- motion resonant orbits with neighboring planets and with every other planet. To 
check orbital stability of such systems, we calculate the crossing time of systems whose 
planets are in 8:7:6 chain resonance, repeatedly (case5). For example, in N = 4 case, a 2nd 
innermost planet and a 3rd one are in 8:7 resonance, and mean-motion ratios of planets are 
28:24:21:18. When planets are in 8:7:6 chain resonance, the pairs of every other planets are 
in 4:3 resonance. Fig. [9] shows that iV cr it ~ 6. This is larger than N CT n = 4 of 7:6 single 
resonance and A^ crit = 3 of 8:7 single resonance. Although chain resonance effect is unclear, 
there is possible that planets would be stabilized by 4:3 every other resonance (iV cr it > 7 in 



10 5 M Q planets) or that it would reduce crossing time dependence on iV like the case of 
mixture of planets in retrograde and prograde orbits (Smith and Lissauer 2009). 

5. Conclusions 

We have investigated the crossing time (t croS s) of resonant systems by numerical 
simulations for 4:3, 6:5, 7:6 and 8:7 resonances. The crossing time of non-resonant plants 
decreases continuously with increasing the total number of planets N. In the case of 
resonant systems, however, while t croS s is comparable to that in non-resonant systems for 
large N, it abruptly changes for N < N cr[t . In that case, the resonant systems are stable 
during the simulation time (10 s Kepler time). We examine 4 cases of different resonances 
with changing the orbital separation in mutual Hill radii K and planetary mass M p , and 
1 case of chain resonances. When K or M p is fixed, iV cr i t increases as p value of p + l:p 
decreases. When planets of mass 1O~ 5 M are in 6:5 resonances, N crit = 8, when planets 
of mass 1O _5 M are in 7:6 resonances, iV crit = 4, when planets of mass 1O _5 M are in 8:7 
resonances, N a n = 3, and when planets of mass 10 _4 M Q are in 4:3 resonances, iV cr i t = 7. 
Observed planets are not always in a chain of the same resonance. We calculate 8:7:6 chain 
resonant case. We find that iV cr it = 6 of 8:7:6 resonance case is larger than either of 8:7 and 
7:6 single cases. 

In this paper, we set the innermost planet at 0.1AU, which is usually too far for tidal 
effect to be particularly important. However, many detected planets reside quite a bit closer 



to their host stars. As can be seen in the numerical simulations of Terquem and Papaloizou 
(2007), Papaloizou and Terquem (2010), and Papaloizou (2011), tidal dissipation affects the 
stability of resonant systems due to two effects. One is eccentricity damping. Another is 
change of amplitudes of resonant libration. According to Papaloizou and Terquem (2010), 
the timescale of eccentricity damping due to the tidal dissipation is 



The parameter Q' = 3Q/2k2, where Q is the tidal dissipation function and &2 is the 
Love number. Since planets are stabilized when e-damping timescale is shorter than the 
crossing time in gas free conditions, tidal dissipation stabilizes planets a < 0.02AU for our 
closely spaced systems. Tidal dissipation often increases amplitude of resonant libration 
in association with an increasing period ratio, Pj+i/Pj (Terquem and Papaloizou 2007; 
Papaloizou 2011). Although tides help eccentricity damping in most cases, planets having 
large resonant amplitudes would easily move away from resonant configurations due to tides 
and cause instability. 

Many planets in observed multi-super Earths systems are orbiting near the central 
star. It is suggested that close-in super-Earths are formed through orbital migration of 
protoplanets and stopped near the disk inner edge due to resonant trapping (Terquem and 
Papaloizou 2007; Ogihara and Ida 2009). Even in 6:5 resonance, our simulation ( §4.21) 
shows the system can hold 8 planets stably. If a planet is in 2:1 mean- motion resonance, the 
orbital separation is larger than 10 Hill radius for a planet M p < 1.4 x 1O _4 M . Observed 
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systems composed of a few planets in 2:1 resonance are stable over 10 s Kepler time since 
these system can hold over 8 planets stably. The same goes for satellite systems such as 
Galilean satellites. On the other hand, observations find many systems that are not in a 
mean-motion resonance (Fabrycky et al. 2012). The systems of super-Earths far from the 
resonances could be formed by the scenario by Ogihara and Ida (2009) that planets more 
than the critical number are once trapped in resonance in gaseous stage and cause orbital 
instability after gas depletion. 



When there is eccentricity damping force, a planetary system is stable (Iwasaki et al. 
2001, 2002), provided that e-damping timescale (tdamp) is shorter than the crossing time 
in gas free case (£* ross ). Since we adopt the exponential decay of gas, tdamp exponentially 
increases and eventually exceeds tcross- We define tdrag as the timescale for tdamp to become 
longer than t* ross . Using the formula by Iwasaki et al. (2002), we find 



A. Timescale of Gas Depletion 



drag 



2.306 (K - K diag (t = 0))r dl 



(Al) 



where 




(A2) 



b and c are defined in Eq. [TJ 
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Table 1: Cases of our calculations. 





p + l:p resonance 


planetary mass M p 


Hill separation K 


critical number A/ crit 


easel 


6:5 


10~ b M Q 


6.450 


8 


case2 


7:6 


1(T 5 M 


5.456 


4 


case3 


8:7 


1(T 5 M Q 


4.727 


3 


case4 


4:3 


1O- 4 M 


4.715 


7 


case5 


7:6 and 8:7 


10~ 5 M„ 


4.727 and 5.456 


6 
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Figure Captions 

Fig. [TJ - An example of resonant trapping is shown. Time evolutions of semi-major 
axes. The total number N is equal to 5. We form exact resonance situations by 5000 
yr simulation of orbital migration before the reduction of the gas density at t = 0. 
Planets migrate inward and are trapped at 7:6 resonance orbits. At t = 0, we start gas 
depletion with depletion timescale Td ep = 10 4 yr. The crossing time is about 73000 yr 
(~ 2.3 x 10 6 T Kcp ). 

Fig. [2J - Time evolutions of resonant angle of 5 planets in 7:6 resonances. In this case, 
the orbital instability causes at about 73000 yr. Planets are labeled from the innermost. 
The upper figure shows resonant angle <pi (cross), y? 2 (square), ip 3 (circle), and y? 4 (triangle). 
The lower figure shows resonant angle ip[ (cross), ip' 2 (square), ip' 3 (circle), and <p' 4 (triangle). 

Fig. [3j - Time evolutions of resonant angles of non-resonant planets. This figure is the 
case of A = 5 and K = 5.456. The orbital instability causes at about 4200 yr. The upper 
figure shows resonant angle </?i (cross), y9 2 (square), y2 3 (circle), and y?4 (triangle). The lower 
figure shows resonant angle ip[ (cross), <p' 2 (square), <p' 3 (circle), and <p' 4 (triangle). 

Fig. ID - Crossing times t cross normalized by Kepler time of the innermost body versus 
the number of planets N. The circles represent the crossing time of 6:5 orbits whose initial 
longitudes are chosen randomly (100 cases for each N). The solid curve is a least square 
exponential fit to the circles, logt croS s = 6.21/A + 3.39. The triangles are the crossing 
time of 6:5 orbits (10 cases for each N) whose initial conditions are generated by orbital 
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integrations including migration (rd cp = 10 4 yr). The square symbols are Td cp = 10 3 yr 
and cross symbols are Td cp = 10 2 yr. Three dotted lines represent tdrag (Eq. IA1I) of the 
outermost planet for Td ep = 10 4 , 10 3 , and 10 2 yr from the top. In the easel, N CTit is equal 
to 8. 

Fig. [5j - Same as Fig. HI 7:6 resonances (case2). The circles represent the crossing 
time of 7:6 orbits whose initial longitudes are chosen randomly (100 cases for each N). The 
solid curve is a least square exponential fit to the circles, logt cross = 5.65/A" + 2.76. The 
triangles are the crossing time of 7:6 orbits (10 cases for each N) whose initial conditions 
are generated by orbital integrations including migration (rd cp = 10 4 yr). In the case2, iV crit 
is equal to 4. The dotted line shows tdrag (Eq. IAip of the outermost planet. 

Fig. [6j - Same as Fig. HJ 8:7 resonances (case3). The circles represent the crossing 
time of non-resonant orbits. The triangles are the crossing time of resonant orbits (3 
runs for each). The dotted line shows idrag of the outermost planet. The solid curve is 
logtcross = 4.91/ N + 2.19. Here, N ciit is equal to 3. 

Fig. [3 - A^crit versus the orbital separation in mutual Hill radii K in the case of 
M p = 10" 5 M o . 

Fig. EJ - Same as Fig. [4j 4:3 resonances (case4). The circles represent the crossing time 
of non-resonant orbits and the solid curve is a least square exponential fit to the circles, 
logtcross = 7.37 /N + 2.05. The triangles are the crossing time of resonant orbits (3 runs for 
each). The dotted line shows tdrag of the outermost planet. In case4, A^t is equal to 7. 



Fig. [9j - Same as Fig. HI but for 8:7:6 resonances (case5). The circles represent the 
crossing time of non-resonant orbits and the solid curve is a least square exponential fit to 
the circles, logt cross = 11.85/-/V + 1.81. The triangles are the crossing time of resonant chain 
orbits (3 runs for each). In case5, iV cr it is equal to 6. 



-32 - 



0.2 
0.18 
0.16 
0.14 
0.12 h 



f = 1 



7:6 
7:6 



7:6 



7:6 



_L 



f g = -exp(t/1 OOOO yrj 



T 



i 



-5000 -4000 -3000 10000 20000 30000 40000 50000 60000 70000 80000 

Time [yr] Time [yr] 



-33 - 



+■ 

AD . 



: 1 1 1 1 | 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 1 1 I I '.'..^ G 




10000 20000 30000 40000 50000 60000 70000 80000 
Time [yr] 



Fig. 



2.- 



-34- 




Time [yr] 

Fig. 3.- 



casel (6:5, K=6.450) 



8 p-f W *|TT 




2 - 



3 5 8 10 20 30 50 
N 

Fig. 4.— 



case2(7:6, K=5.456) 
8 F + 4 I I 1 

7 - 



JS'A" "A" 




3 4 5 6 10 30 50 



Fig. 5- 



case3(8:7, K=4.727) 





CO 


_TII I 


1 1 - 




7 








Q- 
0) 


6 








w 


5 








o 
+ _p 

o 


4 










3 


1 ! 








2 




o 

1 1 T 


1 1 



3 4 5 10 30 50 

N 



Fig. 6.- 



-36 - 




Fig. 8.- 



-37- 



case5(8:7 and 7:6 resonances) 



8 p-fr — I 




3 5 6 78 10 
N 



Fig. 9.- 



